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ABSTRACT 

A three-dimensional Monte Carlo code for modelling radiation transport in Type la 
supernovae is described. In addition to tracking Monte Carlo quanta to follow the 
emission, scattering and deposition of radiative energy, a scheme involving volume- 
based Monte Carlo estimators is used to allow properties of the emergent radiation 
field to be extracted for specific viewing angles in a multi-dimensional structure. This 
eliminates the need to compute spectra or light curves by angular binning of emergent 
quanta. The code is applied to two test problems to illustrate consequences of multi- 
dimensional structure on the modelling of light curves. First, elliptical models are 
used to quantify how large scale asphericity can introduce angular dependence to 
light curves. Secondly, a model which incorporates complex structural inhomogeneity, 
as predicted by modern explosion models, is used to investigate how such structure 
may affect light curve properties. 

Key words: radiative transfer - methods: numerical - supernovae: general 



1 INTRODUCTION 



Despite decades of study, type la supernovae (SNIa) con- 
tinue to be an active topic for astrophysical research. The 
accepted physical mechanism for these events, that they are 
the result of thermonuclear explosions of degenerate ma- 
terial (Hoyle & Fowler 1960) in white dwarf stars, is very 
well established but owing to the considerable complexity of 
this process many aspects of these spectacular events remain 
poorly understood. 

Observationally, it is clear that all SNIa are not the 
same - there is significant diversity in both brightness and 
decay timescale (Phillips 1993). Aside from its direct rele- 
vance to the study of SNIa themselves, understanding this 
diversity has important implications to other branches of 
astrophysics - in particular to cosmology where inferences 
about the properties of distant SNIa play an important role 
in probing the rate of expansion of the Universe. Some corre- 
lations between supernova properties are already well estab- 
lished from observations of nearby SNIa (e.g. Phillips 1993). 
However, it is becoming clear that the observed diversity is 
not adequately described by a single parameter (Benetti et 
al. 2004, 2005) and that unraveling this diversity relies on a 
combination of careful observational study and sophisticated 
theoretical modelling. 

In recent years, there has been rapid development in 
the sophistication of numerical modelling of the hydrody- 
namics of SNIa explosions. In particular, while in earlier 
work supernovae explosions were modelled using first one- 



dimensional (ID) computer codes (e.g. Nomoto, Thielemann 
& Yokoi 1984; Hoflich, Wheeler & Thielemann 1998) and 
later two- dimensions (e.g. Miiller & Arnett 1986; Niemeyer, 
Hillebrandt, Woosley 1996), the most up-to-date simula- 
tions are fully three-dimensional (3D, e.g. Reinecke, Hille- 
brandt & Niemeyer 2002; Gamezo et al. 2003; Ropke 2005; 
Ropke et al. 2006) . These multi-dimensional models are cru- 
cial for the understanding of realistic flame propagation and 
hence nucleosynthesis in SNIa explosions and have clearly 
demonstrated that the earliest ID models greatly under- 
estimate the likely complexity of real SNIa. This naturally 
raises the question of whether and to what extent this multi- 
dimensional complexity may be responsible for the observed 
diversity of explosions. 

To connect hydrodynamical explosion models and ob- 
servations of real SNIa light curves or spectra requires mod- 
elling of radiation transport in the supernova. The majority 
of the light escaping from a supernova explosion at around 
optical maximum originates from energy deposited in the 
ejecta by the absorption and Compton scattering of 7-rays 
emitted by radioactive isotopes. Given that all SNIa explo- 
sion models are significantly optically thick at the earliest 
times (< 1 day), the emission, transport and deposition of 
radiation need only be followed for times after most of the 
complex dynamics have ceased and the ejecta is in near- 
homologous expansion; this greatly simplifies the radiation 
transport calculation and allows it to be decoupled from the 
hydrodynamical simulation. 

The level of sophistication of modern radiation trans- 



2 S. A. Sim 



port computations for supernovae has gradually developed 
to keep pace with the development of the hydrodynami- 
cal explosion models: for earlier applications (e.g. Branch 
et al. 1981, Lucy 1987, Mazzali & Lucy 1993) simple one- 
dimensional calculations were sufficient but there is now 
growing interest in fully three-dimensional simulations. In 
particular, Lucy (2005) has recently described and tested a 
sophisticated Monte Carlo approach to 3D, time-dependent 
radiation transport. This method holds much promise for 
a wide range of applications and is already being adopted 
and developed in several contexts - for example, Kasen, 
Thomas & Nugent (2006) have developed a code based on 
that described by Lucy (2005) which incorporates a treat- 
ment of polarisation and lifts the assumption of grey opacity 
used in the earlier code. Independently, Maeda, Mazzali & 
Nomoto (2006) have also adopted Lucy's techniques and uti- 
lized them in multi-dimensional simulations of light curves 
for core collapse supernovae. 

In this investigation, Monte Carlo calculations using 
methods similar to those of Lucy (2005) are used to investi- 
gate aspects of the influence of the 3D structure of SNIa on 
their observable properties. In Section 2, the computer code 
used for the calculations is described. Although similar to 
the codes described by Lucy (2005) and Kasen et al. (2006), 
there are key differences, specifically in the means by which 
the observable light curves are extracted from the Monte 
Carlo simulation. This code is then used to investigate the 
effect of departures from spherical symmetry for two physi- 
cally distinct cases. First, in Section 3 toy models are used 
to investigate the influence of large scale (low-mode) asym- 
metries on observable light curves. Secondly, a representa- 
tion of a modern 3D explosion model is used to investigate 
the implications of the complex, relatively small-scale inho- 
mogeneity predicted (Section 4). The emphasis here is on 
identifying and interpreting differential effects between 3D 
and ID radiative transfer and thus several simplifying as- 
sumptions are made in the treatment of the micro-physics. 
Conclusions are drawn in Section 5. 



2 METHOD 

The calculations presented in this paper have been per- 
formed using a Monte Carlo code which follows the prop- 
agation of radiative energy in 3D as a function of time. This 
code is closely based on that described by Lucy (2005) and 
thus only a brief summary of the operation of the code is 
given here, with particular emphasis on the the departures 
from the approach described by Lucy (2005). The nomencla^ 
ture used by Lucy for the Monte Carlo quanta (radioactive 
pellets, 7-packets, r-packets) is adopted throughout. 

2.1 Conceptual summary of the code 

To obtain light curves for a pre-specified supernova model, 
the code undertakes the following steps. A computation do- 
main is assigned which is large enough to encompass the 
physical extent of the model. This is divided into a number 
of grid cells each of which is assigned an initial mass-density 
based on the input supernova model. Pellets of radioactive 
material are placed in these cells, again in accordance with 
the chosen model. The radiative decays of these pellets by 



emission of 7-rays and the subsequent propagation and ther- 
malisation of the 7-packets is followed via a Monte Carlo cal- 
culation. A grey-opacity treatment is adopted for the prop- 
agation of radiation (r-packets) in all other wave-bands. 

During the Monte Carlo calculation, estimators are de- 
termined for various properties of the radiation field in every 
grid cell (see below). No information is retained regarding 
the variations of these properties within individual grid cells. 

After the Monte Carlo simulations are complete, the 
estimators arc used to determine observer frame cmissiv- 
ities. A formal solution of the radiative transfer equation 
is then performed using these emissivities to determine the 
emergent radiation field. The code currently produces two 
classes of output; "bolomctric" (ultraviolet-optical- infrared, 
UVOlR)-light curves which are obtained from the behaviour 
of the r-packets, and 7-ray spectra and light curves which are 
derived from the 7-packets. For the applications described in 
Sections 3 and 4, only the r-packet light curves are needed, 
however the analogous means by which 7-ray properties are 
obtained are also described here for future reference. 

2.2 Tracking of quanta and the grid 

Following Lucy (2005), a regular Cartesian grid which ex- 
pands with time is used. However, we allow the grid to ex- 
pand continuously rather than in discrete jumps at the end 
of a each time-step. Since the supernova ejecta is assumed to 
be in homologous expansion, this adds little complexity to 
computing the propagation of the Monte Carlo quanta and 
removes the need to check whether they skip across bound- 
aries due to the modification of the grid at the start of each 
time-step. In the limit of small time-steps, this modification 
has no effect on the results of the calculation. 

2.3 Extraction of spectra and light-curves 

The code described by Lucy (2005) obtained light curves 
by directly counting the number of quanta that escaped the 
computational domain during each time-step. However, as 
discussed by Lucy (2005), that approach is not the most 
efficient with regard to minimizing Monte Carlo noise. Sig- 
nificantly higher quality spectra and light curves can be ob- 
tained by using the paths of the Monte Carlo quanta to com- 
pute estimators for the emissivities in the model which can 
then be used, post hoc, to obtain the intensity via a formal 
solution of the radiative transfer equation. Methods using 
such an approach have previously been used and shown to 
be highly successful (see e.g. Lucy 1999). 

For spherical models - such as the test case used by 
Lucy (2005) - the method of counting packets can always 
be used reliably, its inefficiency being countered by the use 
of a large number of quanta. For complex 3D structures, 
however, this approach becomes unsatisfactory since the an- 
gular dependence of the escaping radiation field can only be 
addressed by angular binning of the quanta. This has fur- 
ther negative impact on the signal-to-noise of the computed 
spectra and rapidly becomes prohibitive if more than a few 
angle-bins are to be considered. 

In contrast, the use of emissivity estimators and a for- 
mal solution of the radiative transfer equation allow the 
spectra and light curves to be computed correctly for in- 
dividual lines-of-sight to the supernova without introducing 
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extra Monte Carlo noise. Thus such methods are strongly 
favoured for the study of any models which depart from 

spherical symmetry. The implementation of these methods ^7™ = '^j^y^ '*t™ ~ 2t'-ftp/c)ds (4) 

for the calculation of UVOIR light-curves and 7-ray spectra 7-paths 
are described in the next two sub-sections. 

and 



2.4 UVOIR light curves 

The light-curve for UVOIR radiation as viewed by a distant 
observer in direction n is given by: 



^obs (^obs 



(1) 



Here, the integral is performed over the plane perpendicular 
to n and Ioc{d,tohs,n) is the emergent intensity of a ray 
which is destined to reach the observer at time tabs- The 
intensity depends on d, the impact-vector of the ray given 
hj d — r X h where r is any position on the ray trajectory. 
Lohs (iobs, n) is the luminosity the observer would imply were 
they to assume that the emission was isotropic. The true 
(unobservable) luminosity of the supernovae, L, is given by 



i(tobs) = — 



all n 



(2) 



This L is the quantity obtained by directly summing all 
emerging packets in a Monte Carlo simulation. 

To evaluate equation 1 for a particular viewing direc- 
tion (n) at a particular time (tobs), a large sample of rays 
(typically 10^) are chosen such that they cover the complete 
range of d required by the size of the supernova ejecta. These 
rays are launched simultaneously from a plane perpendicular 
to (n) behind the supernovae having initially zero intensity. 
The time of launch is chosen such that the ray with d = 
crosses the coordinate origin at time tobs- For each ray, the 
emergent intensity is determined by numerically solving the 
radiative transfer equation along its trajectory: 



dJ(r,t',n) 
d^ 



: r}(r,t' , n) — /t(r , t', n)I{d, t' , n) (3) 



where ds is the element of pathlength along the ray while rj 
and K are the observer frame, direction-, time- and position- 
dependent emissivity and opacity, respectively. Since the tra- 
jectory is a light ray, the position r and the time t' arc re- 
lated by dr/dt' = cn. The emergent intensity, /oo(d, tobs, n), 
is the value of I{d,t' ,"(1) at the point where the ray trajec- 
tory finally leaves the supernova ejecta. 

In all the calculations presented here, the uvoiR opacity 
per unit density in a grid cell is assumed to be constant and 
isotropic in the co-moving frame. Thus the observer frame 
opacity, /t(r, t' , n) can be readily calculated. In contrast, the 
emissivity ri{r,t',n) is not known a priori and is obtained 
from a Monte Carlo simulation. 

In the current version of the code, there are two distinct 
UVOIR emissivity source terms: one due to thermalisation 
of 7-ray packets {ri-,) and one due to scattering of UVOIR 
photon packets (r/r). Following Lucy (2005), estimators for 
these emissivities in the co-moving frame for a particular 
grid cell during a particular time step are obtained from the 
Monte Carlo energy packet trajectories using: 



r-paths 



/t^™'e(l — 2v.'h,p/c)ds 



(5) 



In equation 4, the summation runs over all the trajectories 
of 7-ray packets within the cell (which has volume V) during 
the time-step (which has duration Af); ds is the trajectory 
length and e is the packet energy determined in the observer 
frame. The packet is travelling in direction -hp and the ve- 
locity of the ejecta at the mid-point of the trajectory is v. 
The co-moving frame 7-ray thermalisation opacity (/tSj™*) is 
frequency-dependent and includes contributions from both 
Compton scattering and photoelectric absorption. Note that 
the Compton term accounts only for the rate at which 7- 
rays transfer energy to Compton electrons (which arc as- 
sumed to thermalise in situ). In equation 5, the summation 
is over uvoiR packet trajectories and kJ;™^ is obtained from 
the adopted uvoiR opacity coefficient. 

It has been assumed here that the emissivity is isotropic 
in the co-moving frame and that terms 0{v^ /(?) and smaller 
can be neglected. Note that although the co-moving frame 
emissivity {rf™^) is isotropic, the observer frame emissivity 
(r?) is not, owing to the Doppler terms in the transformation 
between frames. 

The scheme described thus far is valid provided that 
the grid cells are individually optically thin. In practice, this 
condition is violated at eaxly times when the ejecta is com- 
pact and dense. At such times, it becomes unacceptable to 
assign a uniform emissivity to each grid cell. This can be 
overcome, however, by weighting each contribution to the 
emissivities with a factor which accounts for the probabil- 
ity of energy absorbed and re-emitted during the related 
trajectory escaping to infinity. The weighting factor used 
here, which must be applied individually to each contribut- 
ing pathlength included in the sums in equations 4 and 5 
is 



w(r(n),gr) = exp(-T(n) — ^Sr) '^^^^^!^ 5- (6) 

2 or 



where 5t = h.hpKrds and T(n) is the total r-packet op- 
tical depth from the mid point of the trajectory ds to the 
edge of the supernova in the direction n. This form of the 
weighting factor is valid provided that either the total opti- 
cal depth across a grid cell is small or that the contributing 
pathlengths (ds) are all small compared to the typical length 
scale on which the physical properties (e.g. mass density) of 
the model vary. Therefore, in calculations which make use 
of these weighted estimators, the pathlengths that energy 
packets can travel in a single step are not permitted to ex- 
ceed a predetermined maximum step size, dsmax- For each 
time step in the calculation, dSmax is set to either one-tenth 
the width of a grid cell or the distance corresponding to 
5t = 0.1 in the densest grid cell; the larger of these two val- 
ues is chosen. Furthermore, when the weighted emissivities 
are used, the opacity term in equation 3 must be neglected 
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Figure 1. Comparison of UVOIR-Iight curves computed for the 

test model described in Section 2.5. The solid line shows the light 
curve computed by Lucy (2005). The open diamond symbols indi- 
cate the light curve obtained by counting emergent packets with 
the code developed here. The broken line is the UVOIR-Iight curve 
computed using the estimators defined by equations 4 and 5. The 
star symbols show the light curve obtained using the weighted 
estimators (Section 2.4, equation 6). 

since the probability of the emitted energy escaping has al- 
ready been addressed. 

This weighting overcomes the problem of having op- 
tically thick colls at early times but it has two significant 
drawbacks. First, the required computations of r are time 
demanding and lead to a substantial increase in code ex- 
ecution time. Secondly, since the weighting factor is angle 
dependent, separate sets of estimators are required for each 
viewing angle that is to be investigated (in contrast, the un- 
weighted estimators are independent of n). Therefore, this 
method should only be used for computing the light-curve at 
early times - when the density in the grid cells is high - and 
be replaced with the unweighted scheme once the opacity of 
individual cells has dropped sufficiently. 

2.5 Test calculations of the UVOIR light curve 

To test the implementation of the method described above, 
the code has been applied to the test model used by Lucy 
(2005). This is a spherical model, based on that used by 
Pinto & Eastman (2000a), having total mass M = 1.39M0, 
^®Ni mass MseN; = 0.625Mq and a maximum velocity of 
lO'* km s~^. The distribution of Ni is centrally peaked. 
Following Lucy (2005), the grey absorption cross-section of 
0.1 cni'^ g~^ is adopted for uvoiR radiation, and the pho- 
toelectric absorption coefficient is taken from Ambwani & 
Sutherland (1988), adopting a mean value of Z = 14. A 
100'' Cartesian grid is adopted for this test and 2.5 x 10® 
Monte Carlo packets were used in the calculations. 

The light curve obtained by Lucy (2005) for this model 
is shown in Figure 1 (solid line) . This agrees very well with 
the light curve computed here by directly counting the num- 
ber of Monte Carlo packets escaping this model (diamond 
symbols in Figure 1). It also matches well with the light- 
curve computed for a particular viewing angle using the 
weighted estimators described in Section 2.4 (the particu- 



Figure 2. The maximum UVOIR optical depth across a grid cell 
as a function of time in the test model described in Section 2.5. 

lar viewing direction was chosen randomly here). However, 
when the estimators are not weighted (i.e. if equations 4 and 
5 are used directly), the computed light curve is significantly 
overestimated at early times (< 20 days, in this case). As 
discussed in Section 2.4, this is because the grid cells of the 
model are optically thick at early times; the ma^ximum op- 
tical depth across one grid cell is plotted, as a function of 
time, in Figure 2. This shows that significant errors arise if 
the unweighted estimators are used when this optical depth 

is Tcell > 1. 

2.6 7-ray spectra and light curves 

The 7-ray spectrum is obtained following the same principles 
as that used for the UVOiR light curve described above. The 
important difference being that the treatment of 7-rays is 
fully frequency dependent. 

To deal with the frequency dependence, a grid of fre- 
quency points is used to divide the spectrum into small fre- 
quency intervals. One point in this grid is set to the rest 
frequency of each of the radioactive 7-ray emission lines in 
the range of interest and the remaining points are spaced log- 
arithmically between. To compute the spectrum, the same 
scheme of solving the radiative transfer equation along a set 
of rays trajectories through the model is used. Here, how- 
ever, the radiative transfer along each trajectory is com- 
puted multiple times, once for each frequency point in the 
frequency grid, thereby determining the emergent radiation 
field as a function of both frequency and time. 

As for the grey-computations described above, the opac- 
ity term in the radiative transfer equation is known (the 
sum of Compton and photoelectric terms) . There are again 
two emissivity terms which need to be considered. The first, 
direct emission of 7-rays by radioactive decay can also be 
expressed analytically for every grid cell in terms of the 
half-lives of the radioactive isotopes and their initial con- 
centrations in the cell. 

The second emissivity term is due to Compton down- 
scattering. The treatment of this term requires that both the 
angular- and frequency-dependence of the Compton process 
be considered. It is determined via a set Monte Carlo es- 
timators (one per frequency interval per grid cell per time 
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step); the estimator for the frequency interval i in the grid 
cell j for timestep k is given (in the observer frame) by 



Vijk = 



VjAviAtk 



E 



7-paths 



(7) 



where the sum runs over all 7-ray trajectories which lie in 
cell j (having volume Vj) and have frequency appropriate for 
scattering into the frequency interval i (of width Av) during 
the timestep k (which has duration Atk)- Tie is the number 
density of target electrons in the coll, e is the observer- frame 
energy of the 7-ray packet and (da/dfi)^,^^ is the Compton 
differential cross-section for scattering into the direction of 
the line-of-sight (n) in the observer frame. This cross-section 
depends upon the angle between the trajectory and h and 
is dotorminod by applying the Klcin-Nishina formula for the 
cross-section in the co-moving frame. 



3 APPLICATION 1: AN ELLIPSOIDAL MODEL 

In this section, two toy models of elliptical supernovae are 
used to investigate possible observational consequences of 
largo scale asphoricity in supernova explosions. Such an in- 
vestigation is motivated by observational evidence for global 
asphericity in SNIa obtained via polarimetry (see e.g. Howell 
et al. 2001; Wang et al. 2003). The origin of this asphericity 
is not well known: mostly likely it is determined by the de- 
tails of the explosion process itself but may have its roots in 
the properties of a rapidly rotating progenitor. Here, how- 
ever, the objective is not to gain insight to the physical origin 
of such a geometry but rather to study how it might effect 
both the amplitude and shape of observed light curves in 
comparison with spherical explosions. 

Earlier calculations of radiation transport in elliptical 
supernovae have been discussed by Hofiich (1991). In that 
work, a multi-dimensional Monte Carlo code was also used. 
However, the treatment of energy packet generation and 
emission was simplified via the use of a paramctcriscd pho- 
tosphere in contrast to the full treatment of 7-transport and 
deposition employed here. 



3.1 The model 

A simple elliptical model has been constructed, closely re- 
lated to the spherical model used as a test case in Section 2. 
The adopted model has the same total mass and same ^®Ni 
mass as the spherical model. It is also assumed to be in ho- 
mologous expansion and to have uniform density. However, 
the maximum velocity (and hcucc spatial extent) is taken to 
be smaller in the ^-direction than in the x- and 3/-directions 
(symmetry is still Eissumed under rotation about the a-axis) . 
Such a model is intended as a simple description for cases 
in which either the explosion mechanism or the properties 
of the progenitor lead to a large-scale (low angular mode) 
departure from sphericity. 

Two particular realisations of the model will be con- 
sidered here. For both the maximum velocity in the x- and 
j/-directions was kept at 10^ cm s~^. The models differ in the 
chosen maximum velocity in the ^-direction: for one, this ve- 
locity was fixed at 5 x 10* cm s~^, thereby giving an axis 
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time (days) 

Figure 3. Computed light curves for the elliptical model with 
axis ratio 2:1. The solid line sliows the anglc-avcragcd light curve, 
the diamonds show the results for viewing down the long axis and 
the stars for viewing down the short axis. 



ratio of 2:1; while for the other 8 x 10® cm s~^ was adopted 
to give an axis ratio of 5:4. The axis ratio of 2:1 may be 
regarded as extreme - it is comparable to the axis ratio that 
might be present in a very rapidly rotating progenitor (see 
e.g. Yoon & Langer 2005) but there is little evidence to sug- 
gest that this would be preserved during an explosion. As an 
extreme case, however, this model is useful for providing a 
clear indication of the sense in which asphericity can affect 
the light curve. The second ratio adopted (5:4) is compa- 
rable to the Ri 20 per cent asphericity found by Howell et 
al. (2001) for SN1999by and is therefore more likely to be 
indicative of the propertied of real SNIa explosions. 

As discussed in Section 1, the interest here is in prob- 
ing the effect of a 3D treatment of the radiative transfer 
compared to spherically symmetric ID calculations and so 
the simplifying assumption of a constant, grey-uvoiR opac- 
ity is retained. For all the calculations discussed in the sec- 
tion, this grey-uvoiR absorption cross-section remains fixed 
at cr = 0.1 cm^g~^. 



3.2 UVOIR light curves 

Light-curves have been computed for observers viewing the 
ellipsoidal supernovae from infinity along both the major- 
and minor-axes. These were computed using the weighted- 
estimators method described in Section 2.4. The two light 
curves for the model with axis ratio 2:1 are plotted in Fig- 
ure 3, along with the angle- averaged light curve for the 
same model (i.e. the light curve obtained from the arith- 
metic mean of the light curves seen by a large number of 
observers from random viewing angles). 

When viewed down the short-axis, the light curve is 
considerably brighter than as observed down the long-axis. 
Around maximum light, the difference in brightness is ap- 
proximately a factor of 3.5. At later times, as the supernova 
becomes less optically thick, the difference becomes smaller 
approximately a factor of 2.5 at around 50 days. The light 
curve peaks slightly earlier when viewed down the long cixis 
(at ~ 15.3 days compared to ~ 16.4 days if viewed down the 
short axis). Also, the light curve decays more slowly if viewed 
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Figure 4. Computed light curves for the eUiptical model with 
ELxis ratio 5:4. The solid line shows the angle-averaged light curve, 
the diamonds show the results for viewing down the long axis and 
the stars for viewing down the short axis. 

down the long axis; this is characterised by the bolometric 
AMi5-value!4J of 0.55 (long-axis) and 0.64 (short-axis). 

Qualitatively similar, but quantitatively smaller differ- 
ences are seen in the light curves computed for different 
viewing angles using the model with the axis ratio 5:4. For 
this model, representative light curves are show in Figure 4. 

The scale of the angular variation is determined by the 
interplay of several effects. Consider viewing an opaque el- 
lipsoid with axis ratio 2:1 and uniform surface brightness as 
characterised by uniform surface temperature. One would 
expect to find the flux to be twice as large when viewing 
along the short axis compared to the long axis, simply due 
to the increase in apparent surface area. The computed ratio 
exceeds this for two closely related reasons. First, the choice 
of uniform density in the models means that the column 
density to any particular 2:1 ellipsoidal surface within the 
model is always less when viewing down the short axis than 
the long. Thus the optical depth to the surface is smallest 
along the short axis making the light curve brighter when 
viewed from this direction. 

Secondly, contours of constant radiation density (or 
equivalently emissivity) do not exactly follow the 2:1 ellip- 
soidal geometry of the model; in the outer regions there can 
be significant departures, always in the sense that the radi- 
ation energy density is highest at the points of intersection 
by the short axis. This is illustrated in Figure 5 where the 
variation of the r-packet emissivity (ry™') is shown along 
the ellipsoidal axes. A higher radiation energy density on 
the short axis makes the light curve brighter when viewed 
down that axis and so enhances the variation in brightness 
with viewing angle. The origin of this effect also lies in the 
lower column densities along the short axis. At late times, a 
quasi-static description of the radiation field becomes valid 
since packets are not trapped in the model for a significant 
number of timesteps. Under such circumstances, the lower 
opacity along the short axis means that photons preferen- 

^ Commonly used in analyses of light curves following Phillips 
(1993), AMi5 is the change in magnitude between maximum light 
and 15 days after maximum light. 
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Figure 5. Variation of the r-packet emissivity, r?™^, in the el- 
lipsoidal model with axis ration 2:1. j}™^ is plotted along the 
short axis (connected-diamond symbols) and the long axis (his- 
togram). These values were obtained from the Monte Carlo es- 
timators (equation 5) in the grid cells lying closest to the axes. 
For ease of comparison, they are plotted against velocity in units 
of the maximum speed (umax) along the relevant axis; thus if 
j^cmf ^gfg iQ follow the 2:1 ellipsoidal geometry of the model, the 
diamond symbols and histogram would lie on the same curve. 
Results are shown for two different times: the lightly drawn lines 
are for t ~ 12 days while the heavy lines are for t ~ 40 days. The 
t ^ AO day values have been offset upwards by +2 dex. 

tially diffuse in that direction, making the energy density 
on an 2:1 ellipsoidal surface peak at the points of inter- 
section with the short axis. At earlier times, this effect is 
enhanced by the time-dependent nature of the calculations; 
fewer packets manage to reach outer ellipsoidal shells than 
would be predicted in a quasi-static description because they 
have had insufficient time to diffuse far enough. The differ- 
ence in diffusion time means that this affects the energy 
density along the long axis more significantly than along 
the short axis and thus acts to enhance the effect expected 
from the quasi-static case. As time passes, the shapes of the 
contours of radiation energy density evolve slightly and the 
outward decline becomes somewhat less steep (this can be 
seen by comparing the t 12 days and t ^ 40 days results in 
Figure 5). However, throughout the time range considered 
here, a departure from 2:1 ellipticity remains. 

Since both the effects described above are the result of 
angular variations in the optical depth to 2:1 ellipsoidal sur- 
faces, both persist while the ejecta remains optically thick. 
They slowly decrease in strength during the decay phase as 
expansion causes the optical depths to drop; at very late 
times the entire ejecta will become optically thin to UVOIR 
radiation such that the light curve will become independent 
of viewing angle. For the model adopted here, however, this 
nebular phase will not begin until several hundred days later 
than the times considered here. 

The viewing-angle dependency of the light curves com- 
puted from these simple models may have some interesting 
ramifications for understanding the observed properties of 
SNIa light curves. As pointed out by Wang et al. (2003) 
in the context of SN2001el, directional dependence of the 
luminosity as predicted for elliptical models of supernovae 
would lead to dispersion in the observed peak magnitudes 
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(based on the earlier work by Hoflich 1991 and their imphed 
asphericity of « 10 per cent for SN2001el, they speculate 
that this dispersion would be around 0.1 mag). The results 
obtained here support this argument and indicate that if 
the degree of asphericity were larger in some cases (e.g. SN 
1999by; Howell et al. 2001) the spread in the peak magni- 
tude could be greater, w 0.4 mag. 

Furthermore, the full light curves computed here al- 
low this dispersion relative to the known relationship be- 
tween light curve shape and peak luminosity to bo exam- 
ined. This trend, the so-called "Phillips relation" following 
Phillips (1993), expresses the negative correlation between 
peak brightness and the AAfi5-pararnctcr measured from 
observed light curves of SNIa. For both the models consid- 
ered here, the variation of AM15 with viewing angle is in 
the opposite sense to the standard relationship. This is il- 
lustrated in Figure 6 where the six light curves shown in Fig- 
ures 3 and 4 are represented as points in the AA/15-A/peak 
plane. The gradient of the standard Phillips relation (de- 
scribing the mean observed relationship between AAfis and 
Afpeak) in the B-band is plotted for comparison in the Fig- 
ure. 

This effect would lead to a detectable scatter about the 
Phillips relation and thus may have a major role to play in 
understanding the diversity of supernova observations: the 
results plotted in Figure 6 would suggest that if SNIa explo- 
sions were moderately elliptical (such that the model with 
axis ration of 5:4 were approximately applicable), viewing 
angle effects could explain a scatter of several tenths of a 
magnitude about the mean relationship. Significant caution 
must be applied in interpreting this result since the grey- 
treatment adopted here does not allow band-limited light 
curves to be studied for direct comparison with observa- 
tions - quantitative differences may occur if the frequency- 
dependence of the opacity were taken into account. Further- 
more, the models used here have predicted an angular varia- 
tion of the radiation energy density and thus, by implication, 
the temperature of the ojocta. Such a variation further con- 
tradicts the used of a uniform opacity and highlights the 
need for the consideration of more detailed micro-physics. 
Also, the particularly simple model chosen (uniform den- 
sity with centrally concentrated ^^Ni and time-independent 
mean opacity) produces significantly smaller absolute val- 
ues of AM15 than are typically observed - thus further 
work using more realistic models of aspherical supernovae 
are needed. 
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APPLICATION 2: 
MODEL 



AN INHOMOGENEOUS 



Modern SNIa explosion models predict complex, three- 
dimensional sub-structure (e.g. Reineckeet al. 2002; Gamezo 
et al. 2003; Ropke 2005; Ropke et al. 2006) within the ex- 
plosion. In contrast, most models used to compute light 
curves for comparison with observation have a smooth, one- 
dimensional density- /composition-profile. 

In this section, the effect of the predicted inhomogene- 
ity oil model light curves will be investigated. In principle, 
there are two classes of effect which are of interest here: first, 
inhomogeneity can change the viewing-angle averaged light 
curve compared to that of a spherically symmetric model. 




Figure 6. The figure shows the peak bolometrlc magnitude 
Mpeak versus AM15 for the light curves computed from the el- 
liptical SNIa models. The two points marked with diamonds in- 
dicate the light curves obtained for viewing the model with axis 
ratio 2:1 along its short and long axes. The stars indicate the 
light curves corresponding to the same two viewing angles for 
the model with axis ration 5:4. The triangles are for the angle- 
averaged light curves of the two models. The solid line indicates 
the gradient of the observed relationship between AM15 and the 
peak magnitude as measured in the B-band (Phillips 1993). The 
normalisation of the observed relationship has been adjusted to 
approximately match the computed angle-averaged values. 



Secondly, inhomogeneity could also lead to angular depen- 
dence of the light curve - the degree of angular dependence 
may be resisonably expected to depend on a combination of 
the dynamic range of the inhomogeneity in density and the 
physical length scale of the variations. The length scale is 
relevant since if only small scale inhomogeneity is present 
it will tend to be averaged out when integrating over solid 
angle and so not introduce significant angular variation of 
the light curve. 

Given that the code used here adopts several simplifying 
assumption - most importantly, perhaps, that of a grey ab- 
sorption coefficient for UVOIR radiation - the emphasis here 
is on understanding and assessing the differential effect of in- 
troducing 3D structure by comparing with an equivalent ID 
model using a fixed set of well-understood approximations. 
Such a calculation is an important first step in understand- 
ing the role of 3D structure and is a useful starting point 
for further work where the micro-physics is improved (see 
Section 5 for further discussion). 



4.1 The model 

The model used for this test calculation is based on the 3D 
explosion model computed by Ropke (2005). For this model, 
Ropke (2005) followed the hydrodynamics of an exploding 
white dwarf star (with total mass 1.4 Mq) for 10 seconds 
on a schematic grid of 512^ Cartesian cells. Only one spatial 
octant was simulated and symmetry under reflection was 
assumed to describe the remaining octants (thus numerically 
only 256^ grid cells were used). The distribution of mass- 
density and mass- fraction of iron-group elements in each grid 
cell at the end of their simulation was made available for this 
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work. Light curves have previously been simulated from ID 
representations of this model by Blinnikov et al. (2006). 

In order to make the calculations hero tractable - in 
terms of both computer memory required and photon statis- 
tics in each grid coll - the model adopted here uses only 170^ 
grid cells. It was obtained from the 512^ grid by first remov- 
ing the outermost cell from both ends of the grid in each of 
the three Cartesian directions resulting in a 510^ grid. The 
mass in the cells removed in this process was negligible. The 
resolution was then reduced by a factor of 3 by subdividing 
the 510^ grid into 3^ blocks and replacing each block with 
a single cell whose density was equal to the mean density of 
the original 27 cells. 

Given the grey treatment of UVOIR radiation currently 
adopted in the code, it is not necessary to specify the de- 
tailed composition of the material in each grid cell. How- 
ever, it is necessary to specify the initial distribution of ^®Ni 
which provides the source of radiative energy. The hydrody- 
namics code used by Ropke (2005) provides as estimate of 
the fractional mass of iron-group olcnionts in each grid cell 
but does not give a reliable estimate of the breakdown of 
this material into specific isotopes and elements. The ^®Ni 
mass-fractions used here were obtained by adopting a con- 
stant ratio for the mass of ^^Ni to the total mass of iron- 
group elements in all grid cells. This ratio was fixed to yield 
a total ^®Ni-mass of 0.28 Mq, as derived by Kozma et al. 
(2005) for this model; note that the nucleosynthesis calcu- 
lations described by Kozma et al. (2005) do not produce 
compositional information in sufficient detail to reconstruct 
the full 3D-distribution of ^®Ni-mass owing to the modest 
(27^) number of tracer particles they used in comparison to 
the number of grid cells in the 3D model. 

In order that the differential efi^cct of the 3D structure 
could be assessed, a ID comparison model was made by 
averaging the 3D model over spherical shells, taking care 
to conserve the radial distributions of total mass and -"'^Ni- 
mass. In order to illustrate the degree of inliomogencity in 
the 3D model, Figure 7 shows the velocity distribution of 
density in the spherically averaged model with points indi- 
cating the actual densities of individual grid cells in the 3D 
model. This shows that across most of the velocity range 
of the model, there is a spread in density of at least a fac- 
tor of three between different cells with similar velocities. 
This is comparable to the dynamic range of density for a 
given velocity in the original hydrodynamical model; thus 
one may be confident that the model used here contains a 
fairly reliable representation of the inhomogeneity implied 
by the hydrodynamics. 

4.2 Treatment of opacity 

We continue to adopt a grey UVOiR-absorption cross-section. 
However, to correctly evaluate the influence of 3D structure 

on the light curve, it is necessary to consider that composi- 
tional inhomogeneity would cause the cross-section per gram 
to be a function of position. It goes beyond the scope of this 
paper to undertake full calculations of the composition de- 
pendence of the opacity; instead, a simple one-parameter 
description of the opacity is adopted, following Mazzali & 
Podsiadolwski (2006). They consider the opacity to be de- 
termined by the mass fraction of iron-group elements on 
the assumption that the opacity per gram is a factor of ten 



higher for the iron-group than for lighter elements. Adopt- 
ing the same assumption, the uvoiR-absorption cross-section 
per gram used in this section is given by 

a = 7V(0.9XFe.g,p + 0.1) [cm^ g-^] (8) 

where Xpc-grp is the mass fraction of iron-group elements, 
which varies from cell to cell, and the normalisation factor N 
is chosen such that the mean value of a in the ejecta is fixed 
to < (T >= 0.1 cm^ g~^. Although crude, this parameterisa- 
tion captures the essential physics that the heavy elements 
dominate the mean opacity and only requires the composi- 
tional information which is directly available from the explo- 
sion model (namely, the total iron-group mass). Note that, 
for ease of comparison with the constant-cr calculations used 
in the earlier sections of this paper, the time-dependence 
imposed on a by Mazzali & Podsiadolwski (2006) is not in- 
cluded here. 

In view of the simplifications used - in terms of both 
numerical resolution and particularly the very simple treat- 
ment of opacity - the results obtained here should not be 
regarded as a definitive prediction of the radiation proper- 
ties of the Ropke et al. (2005) explosion model. Rather, the 
emphasis is on using a model which has a reliable representa- 
tion of the degree of inhomogeneity in a real explosion model 
to understand the role played by the complex structure in 
the radiation transport. 



4.3 UVOIR light curves 

The light curve obtained with the 3D model described above 
is shown as the solid histogram in Figure 8; this is the light 
curve averaged over viewing angles. At peak, this light curve 
is dimmer than those plotted in Figures 1, 3 and 4; a con- 
sequence of the lower '"''Ni-mass in the present model. The 
light curve peaks earlier (at t ~ 13 days) as a result of differ- 
ences in the distribution of ''^Ni with velocity: in the model 
used here, the ^^Ni-distribution is less centrally concentrated 
than that adopted in Sections 2.5 and 3. 

To establish the influence of the inhomogeneity on the 
light curve, the ID model obtained by spherically averag- 
ing the 3D model (see above) was also used to compute a 
light curve. This light curve is shown by the dashed curve 
in Figure 8. 

From the Figure, it is apparent that the inhomogeneity 
causes the light curve to be a little brighter at early times 
and very slightly fainter during the decay phase. This oc- 
curs because in the complex 3D medium a fraction of the 
r-packets find trajectories that preferentially pass through 
lower density cells and therefore encounter lower opacities 
than is possible in a spherically symmetric model with the 
same total mass. This allows some r-packets to escape more 
quickly making the light curve brighter during the pre- 
maximum rise phase. It also means that slightly fewer pack- 
ets remain trapped in the ejecta which leads to the marginal 
dimming in the post-maximum phase. The scale of the ef- 
fect in this model is rather modest: at most it amounts to 
a difference of « 0.1 mag and this occurs only at very early 
times, 10 days before maximum light. At around maxi- 
mum light, the effect is limited to a change in magnitude of 
around AM « 0.05 mag (5 per cent in luminosity). 
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Figure 7. The solid line shows the density distribution with velocity of the spherically averaged ID model. The grey points eax;h indicate 
the density and mid-point velocity of a grid cell in the 3D model. All densities are shows for time t = 10 sec. 



In accordance with "Arnett's Rule" (Arnett 1982; Ar- 
nett, Branch & Wheeler 1985) - which states that the emit- 
ted bolometric luminosity is roughly equal to the instanta- 
neous rate of generation of radioactive luminosity at maxi- 
mum light - the increase in peak luminosity resulting from 
the 3D structure is accompanied by a faster rise to max- 
imum; the peak luminosity occurs approximately one day 
earlier in the 3D calculations than the ID case. 

The 3D model was also used to look for viewing-angle 
dependency of the light curve. Light curves were extracted 
for a set of eight randomly chosen viewing angles and 
compared to the angle-averaged light curve. No significant 
(> 1 — 2 per cent) departures from the average light curve 
were found; this can be understood since large scale depar- 
tures from sphericity are not present in this model, partially 
as a result of the assumed reflection symmetries (see above) . 

It is interesting to note that the results from the model 
considered here show that inliomogcncity would act in the 
sense of pushing model results closer to observation; for ex- 
ample, a recent comparison of light curves computed from 
ID models with observations (Blinnikov et al. 2006) does 
suggest that the ID models produce light curves which rise 
too slowly and underestimate peak luminosities (sec their 
figures 8 - 13). Direct comparisons with the calculations of 
Blinnikov et al. (2006) are not possible here owing to the 
considerably more complex treatment of opacity that they 
adopt, but it may reasonably be speculated that their rise 
times would also be shorter if multi-dimensional effects were 
incorporated. 

The scale of the effect here is rather too small to directly 
have major implications for the confrontation of models and 
observations. However, the effect is not completely negligible 
and it is plausible that the calculations presented here may 
underestimate the true scale of this phenomenon; in partic- 
ular, both the reduction in resolution from that used in the 
explosion model and the grossly simplified onc-paramctcr 
treatment of the opacity may both suppress 3D density and 
composition effects that would be present in a complete 
treatment. Furthermore, only one explosion model has been 



-18.4 
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Figure 8. The light curve obtained from the 3D inhomogeneous 
model described in Section 4.1 (solid histogram) and the com- 
parison spherically averaged model (dotted line). The light curve 
computed from the 3D model does not vary significantly with 
viewing angle. 

considered here and this model is constructed from a simu- 
lation describing only one octant (see above). In the future, 
it will be necessary to examine a range of models including 
fully 3D models where the greater complexity may enhance 
the effects. 



5 DISCUSSION 

We have described and tested a new 3D, timo-dcpondont 
Monte Carlo code for modelling radiation transport in SNIa. 
The code adopts the methods presented by Lucy (2005) and 
also incorporates a scheme which uses Monte Carlo radiation 
field estimators to allow obsorvablcs to be extracted for spe- 
cific viewing angles; this approach helps to suppress Monte 
Carlo noise and eliminates the need to obtain intensities by 
angular binning of emergent Monte Carlo packets. 
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This code has been used to investigate two classes of 
three dimensional effects in SNIa models. First, two ellipti- 
cal SNIa toy models were used to investigate how large scale 
asphericity might alter observable light curves. As expected 
from previous simplified treatments (e.g. Hoflich 1991), it 
was found that light curves were brighter when viewed down 
the minor-axis than the major-axis. The brightness enhance- 
ment is largest at eaxly times and becomes smaller during 
the decline phase as the supernova ejecta becomes less opti- 
cally thick. For a model with axis ratio comparable to that 
suggested by polarisation data for real SNIa (Howell et al. 
2001; Wang et al. 2003), the differences in peak brightness 
and light curve shapes between viewing angles is detectable 
and in principle could play a role in the observed scatter 
about the mean relationship between light curve brightness 
and width. 

Secondly, a model with structure based on the results 
of recent three-dimensional hydrodynamical simulations was 
used to study the effect of inhomogeneities in both den- 
sity and composition on light curves. It was found that 3D 
structure could lead to light curves which are both brighter 
at early times and which peaked sooner after the explosion 
that is found from ID models. For the particular model con- 
sidered, the effect was rather modest (« 5 per cent around 
maximum light) but study of other explosion models is re- 
quired to quantify the possible diversity in this effect. 

Considerable further work is required to fully under- 
stand the role of 3D effects on radiation transport in SNIa. 
The use of the grey approximation in the treatment of UVOIR 
radiation is currently the greatest limiting factor to the prac- 
tical applicability of these results. In calculations involving 
a uniform grey absorption coefficient, the line-of-sight opac- 
ity depends solely on the total column density. Under re- 
alistic conditions for SNIa ejecta, the opacity is dominated 
by spectral lines (see e.g. Pinto & Eastman 2000b) and is 
thus a strong function of frequency, velocity gradient, com- 
position and the ionization/cxcitation state of the plasma. 
Indeed, when strong lines dominate, the opacity has little di- 
rect sensitivity to the column density and is instead mostly 
a function of the density of spectral lines in frequency space. 
Furthermore, photon escape from highly opaque material is 
facilitated by frequency redistribution to those regions of the 
spectrum where there are relatively few spectral lines (see 
e.g. Pinto & Eastman 2000b for a detailed discussion); in 
the grey calculations this means of escape is not available 
and all the energy packets must burrow through the imposed 
opacity. Given this, it is likely that the grey approximation 
overestimates the role of geometry in determining photon 
propagation. However, further calculations are required to 
determine and understand this in detail. In the context of 
the inhomogeneous model discussion in Section 4, a more 
realistic treatment of the role played by compositional in- 
homogeneity would be of particular interest; only a crude 
attempt at describing the variation of opacity with compo- 
sition has been used here (equation 8) and a sophisticated 
treatment involving the very real differences between the 
atomic properties of different elements may lead to interest- 
ing effects. 

Fortunately, as discussed by Lucy (2005), the Monte 
Carlo method is well suited to incorporating realistic physics 
in a complex geometry which will allow the interplay of more 
detailed micro-physics and 3D structure to be studied in the 



near future; this would involve a more sophisticated treat- 
ment of the opacity, perhaps similar to the methods em- 
ployed by Bhnnikov et al. (2006) or Kasen et al. (2006). 
Prom such a calculations it would be possible to extract 
not only a bolomctric light curve, but also light curves for 
specific photometric bands, an important step in the pro- 
gression towards confrontation of theory and observation. 

Light curve calculations also need to be performed for 
a wider range of explosion models than the exploratory set 
discussed here. This is necessary in order to quantify how 
the various effects might differ in scale depending on vari- 
ous properties of the explosion and so help to understand 
what role they might play in establishing both the average 
properties and diversity of SN observations. 

Eventually, it would also be valuable to undertake a 
similar assessment of the role of 3D effects in full NLTE 
modelling of UVOiR-spectra; much of the greatest diagnostic 
power lies in spectral modelling including both intensity 
and polarimetry - and thus it is important that commonly 
used modelling assumptions, such as that of spherical sym- 
metry, are studied in the context of modern explosion mod- 
els. 
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